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Abstract 

We use the graded eigenvalue method, a variant of the supersymmetry technique, 
to compute the universal spectral correlations of the QCD Dirac operator in the 
presence of massive dynamical quarks. The calculation is done for the chiral Gaus- 
sian unitary ensemble of random matrix theory with an arbitrary Hermitian matrix 
added to the Dirac matrix. This case is of interest for schematic models of QCD at 
finite temperature. 
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1 Introduction 

By now, it has been firmly established that the spectrum of the QCD Dirac 
operator possesses a number of universal features which can be described by 
chiral random matrix theory (RMT) [1,2]. In particular, the RMT predictions 
agree very well with data from lattice gauge simulations, both for the eigen- 
value correlations in the bulk of the spectrum on the scale of the mean level 
spacing [3-5] and for the distribution and correlations of the low-lying eigen- 
values [6-9]. It is the key assumption of chiral RMT that the matrix elements 
of the Dirac operator in an appropriate energy basis behave as random num- 
bers. This concept has already been very successful in many other areas of 
physics, see the detailed review in Ref. [10]. Thus, it is fair to say that RMT 
approaches can be viewed as thermodynamics for spectral fluctuations and 
related properties. In the context of QCD, deviations from pure RMT statis- 
tics have also been found in the microscopic [11,12] and in the bulk region [5]. 
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These findings provide evidence for the conjecture that lattice QCD may have 
much in common with disordered systems [13,14,11]. In the present work, how- 
ever, we focus on some formal and theoretical aspects of chiral RMT. Hence, 
we shall not discuss in any detail the physical applications but refer to the 
existing literature [15,16]. 

Our aim here is to compute the universal spectral correlations of the QCD 
Dirac operator in the presence of massive dynamical quarks for the chiral 
Gaussian unitary ensemble of RMT. Such a calculation has been done previ- 
ously at zero temperature using orthogonal polynomials [17,18] and the finite 
volume partition function [19,20]. In Ref. [21], the connection between these 
two results was established in the framework of partially quenched chiral per- 
turbation theory. In this paper, we employ the graded eigenvalue method 
[22-24], which is a special variant of the supersymmetry method [25,26], for 
the following two reasons. First, it allows us to extent the calculation to the 
case where an arbitrary deterministic Hermitian matrix is added to the Dirac 
matrix. This is relevant for schematic random-matrix models of QCD at finite 
temperature [27,28]. The standard orthogonal-polynomial method cannot eas- 
ily be applied in this case, because a certain rotation invariance in the space 
of the random matrices is lost. Second, the present problem leads to an in- 
teresting extension of the graded eigenvalue method in the context of chiral 
RMT. In ordinary RMT, this method was developed to calculate spectral cor- 
relations in crossover transitions from regularity to chaos [23,24]. This method 
was then extended to chiral RMT and used to compute the universal spectral 
correlations of the Dirac operator in the quenched approximation, i.e., without 
dynamical quarks [29,30]. For an application of the supersymmetry method 
to the case of one flavor at zero temperature, see Ref. [31]. 

In Refs. [17,32] it was shown that the RMT results are invariant under de- 
formations of the distribution of the random matrix. While we suspect this 
statement to hold also for the results computed in the present paper, a rig- 
orous proof would require an extension of the work of Ref. [33] . We shall not 
address this issue here. 

The outline of this paper is as follows. In Sec. 2, we define the problem and 
outline the main idea for the solution. The graded eigenvalue method is applied 
in Sec. 3. The special case of energies and masses on the scale of the mean 
level spacing near zero, which is also the most interesting case for physical 
applications, is considered in Sec. 4. We conclude with a summary in Sec. 5. 
Technical details of the calculation are discussed in two appendixes. 
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2 Setup of the calculation 



Since this paper is a natural extension of Ref. [29], we shall attempt to use 
a similar notation. The QCD Dirac operator in Euclidean space is defined by 
D = , ~f^d f j i + ig'-f^A^, where g is the coupling constant and the are the gauge 
fields. Note that D is anti-Hermitian. In a random matrix model in a chiral 
basis, the matrix A representing the Dirac operator has the form [2] 



D 



iA 



W + Y 
+ Y 



where W is a square random matrix of dimension N and Y is an arbitrary 
Hermitian matrix. Expression (1) is a schematic model for the QCD Dirac 
operator at finite temperature. The matrix Y represents the effects of the 
temperature on the Dirac spectrum. Its specific form depends on the choice 
of basis states [27,28]. Since we consider an arbitrary Hermitian matrix Y, we 
cover all possible choices of basis states. 



One could also consider the more general problem of a rectangular matrix W, 
giving rise to exact zero modes of the Dirac operator whose number can be 
identified with the topological charge. At Y = 0, this is not necessary since this 
problem is equivalent to introducing additional massless flavors [34] . Although 
it is not obvious, we expect this equivalence to hold also for nonzero Y. We 
hope to address this problem in future work. 

We will be interested in the correlations of the eigenvalues of the matrix A. 
In this paper, we study the chiral Gaussian unitary ensemble (chGUE) appro- 
priate for QCD with three or more colors for which W is a complex matrix 
without any symmetries [34]. The probability distribution of W is given by 

P(W) = - P (W) 1] det(im / - A) , (2) 
yv /=i 

with 

P (W) = exp (-AS 2 tr WW f ) . (3) 

Here, Ay is the number of quark flavors with masses rrif, J\f is a normalization 
factor (see below), and X is a real parameter which will turn out to be equal 
to the chiral condensate at zero temperature. Note that M depends on the 
quark masses. To be precise, the argument of the determinant in (2) should 
have been my + D — nif + iA. For convenience, we have pulled out a factor 
of —% and absorbed it in M . 
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We are interested in the /c-point spectral correlation functions, defined as 
the probability of finding energies in infinitesimal intervals around the points 
x±, . . . , Xk, regardless of labeling. Apart from some trivial contributions involv- 
ing 5-functions [10] they are given by 



Rk(xi, ...,x k )= (--)" I d[W]P(W) n Imtr - J— , (4) 

\ 7T/ J p=1 X p A 



where with e positive infinitesimal. The measure d[W] is simply 

the product of the differentials of the real and imaginary parts of the ele- 
ments of W, i.e., of all independent variables. The integrations extend from 
— oo to +00. The normalization factor M is determined by the requirement 
/ d[W]P(W) = 1. Since both the distribution (2) and the measure d[W] are in- 
variant under unitary transformations of W, only the relative unitary rotation 
between W and Y matters and we can, without loss of generality, write the 
Hermitian matrix Y in diagonal form, Y = diag(yi, . . . ,Un)- Advantageously, 
the /c-point functions can be obtained as 



/ 1 \ k d k 



(5) 



with a generating function given by 

Z k {J) = J d[W]P(W) n det(x p - A) Im l _ - _ - , (6) 

where J stands for J 1; . . . , J^. The starting point of the graded eigenvalue 
method is to rewrite the determinants in (6) as Gaussian integrals over com- 
muting and anti-commuting variables. Note that the distribution P(W) con- 
tains Nf determinants in the numerator. To retain the determinant structure 
of the problem, it is highly desirable to have an equal number of determinants 
in numerator and denominator so that the bosonic and fermionic blocks in the 
super symmetric representation of the generating function have the same size. 
Therefore, we introduce Nf additional determinants in the denominator and 
write, in the large- A" limit, 

Z k (J) = hm Z k {J) (7) 

with 

Udet( imf -A) Im -i4) , (8) 



X 
/=1 
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where the a/ (/ = 1, . . . , Nf) are dummy real variables. The modified normal- 
ization Af depends on the aj and is given by 

Af = J d[W]P (W)J[det(im f -A) Im^-^— -. (9) 

The fact that Eq. (7) holds in the limit N — > oo is proved in App. A. The 
introduction of the dummy determinants in (8) is the main idea of the present 
calculation. As we shall see below, it allows us to use the results of Ref. [29] 
so that the generalization from the quenched approximation to the case with 
Nf > can be obtained with moderate effort. 

We observe that (9) is essentially a special case of (8) with k = 0. To simplify 
the notation, we will compute the generic quantity 

Gj (t) = J d[W\ P (W) n det(f i2 - A) Im _ A) , (10) 

where 7 is a nonnegative integer and t = diag(tn, . . . , t 7 i, in, . . . , t^) is a 
diagonal graded matrix of dimension 27. Both (8) and (9) can be obtained 
from (10) by choosing 7 and t appropriately. 



3 Supersymmetric representation and graded eigenvalue method 



Since we can employ the results of Ref. [29] to compute the function (10) 
very efficiently, we only review the major steps in the derivation. First, the 
determinants in (10) are rewritten as Gaussian integrals over commuting and 
anti-commuting variables which are arranged in a graded (or super) vector ip. 
Then, the integration over W can be performed, resulting in fourth-order terms 
in the -^-variables. These terms can be removed by a Hubbard-Stratonovitch 
transformation at the expense of introducing additional integration variables 
which can be arranged in a complex graded (or super) matrix o. The order of 
the integrations over o and ip can then be interchanged, provided that one is 
only interested in the imaginary parts in (10). This point has been discussed 
in Refs. [29,30,35]. The -^-integration can then be performed trivially. The 
graded matrix a can be written in spherical coordinates as usv with graded 
(or super) unitary matrices u and v and radial coordinates s. The integration 
over u and v can be performed using the supersymmetric generalization of 
the Berezin-Karpelevich integral [36]. It can be viewed as the extension of the 
supersymmetric Itzykson-Zuber integral [22] to complex graded matrices. The 
final result for the function (10) then becomes (see Eqs. (33) through (36) of 
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Ref. [29]) 



C(t)-( 71 ^ ^xpl-iV^trgt 2 ) 



with 



C N {x 1 ,x 2 ) = (2iV£ 2 ) 2 / / d Sl ds 2 2 glS2 2 exp (-iV£ 2 (s 2 + s 2 

J J S~\ ~\~ So 

12 

^ u 2 + s 2 

x J (2iVS 2 Sia ; 1 ) J (2iV£ 2 S2 :r 2 ) Im J] 2 " , > ( 12 ) 

n=l Vn - \ s l ) 

where J and I denote the Bessel and modified Bessel function, respectively. 
In Eq. (11), the symbol trg denotes the graded trace, and 

B,(t 2 ) = n 7 ffi ^ffl with A 7 (x)=f[(x i -x i ). (13) 

Equations (8) and (9) can now be obtained by choosing 7 = k + Nf and 
7 = Nf, respectively, and substituting appropriate values for the entries of t. 
Performing the differentiations according to Eq. (5) (with Z k (J) replaced by 
Z k (J), see Eq. (7)) then yields the Appoint functions. We obtain after some 
algebra 

D / \ f 2 \ k TT 1 1- det[C N (z p ,(q)]p,q=l,...,k+N f , 

Rkix!, . . . , x k ) = - M [ x p hm ; — ^ (14) 



with 





for p = 


1, . . . , k, 




&p—k 


for p = 


k + 1,... 


,k + N f , 


Xq 


for q = 


1, . . . , k, 




im q _ k 


for q = 


k + 1,... 


,k + N f . 



(16) 



Note that Eq. (14), just like Eq. (7), holds in the limit N — > 00 (which is 
the interesting limit for physical applications) but not for finite N. Instead 
of introducing dummy determinants in Eq. (8), there is an alternative way 
to proceed which is exact for finite N. We briefly explain the idea for read- 
ers familiar with the graded eigenvalue method. Without the introduction of 
the dummy determinants in Eq. (8), the transformation from ordinary space 
to superspace leads to a graded (or super) matrix a whose boson-boson and 
fermion-fermion blocks have dimension k and k + Nf, respectively. The trans- 
formation of a to spherical coordinates involves a Berezinian which cannot be 
written as a determinant. The idea now is to enlarge the boson-boson block 
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of a to dimension k + Nf by introducing dummy integration variables in su- 
perspace. Then, the Berezinian resulting from the enlarged u-matrix can be 
written as a determinant which is the prerequisite for expressing the Appoint 
function in form of a determinant. However, we will not discuss this alternative 
way since we are only interested in the limit N — > oo for which the method of 
Eqs. (7) through (9) appears to be more economic. 

Our conventions are such that the support of the spectral density is of order 
(9(1) and the typical level spacing is of order 0(1/N). While Eq. (14) holds for 
all values of the x p and m/, we are particularly interested in the microscopic 
region where the x p and rrif are of order 0(1/ N). We now turn to this limit. 



4 Microscopic limit 



If Xi and x 2 in Eq. (12) are of order 0(1/ N), the integrals can be performed 
in saddle-point approximation in the large- N limit. This was done in Ref. [29], 
and we obtain the Bessel kernel (see Eq. (63) of [29]) 

AT _ x 1 J 1 (2NEx 1 )J (2NEx 2 ) - x 2 J (2NEx 1 )J 1 (2NEx 2 ) 

UjV(£l, X 2 ) — 7ll\ZL g 2 ' 

X l ~ X 2 

(17) 

where E is the only real and positive solution of 

1 N 1 
1= ivS(&/„) 2 + (H/X) 2 (18) 

or zero if no such solution exists [29]. As we shall show below, E = E(Y) can 
be identified with the chiral condensate in the presence of the arbitrary offset 
Y . We now rescale the energies, the masses, and the dummy variables by 2NE 
and define u p = 2NEx p , {if = 2NEmf, and otf = 2NEaf. We thus obtain 
from (14) 

D I \ fOAT^\k I A IT det[C(z p , (q)]p,q=l,...,k+N f , . 

R k (x!,..., x k ) = (2Nz) \[[up\ hm — 19 

\pA ) det[C(a f ,tfig)] f:g=1: ... >Nf 

with a kernel given by 

r, ( \ u iJi( u i) J o( u 2) - u 2 J (u 1 )J 1 (u 2 ) 

l^[U 1 ,U 2 ) - ~2 — 2 , [2U) 

u l a 2 
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where we have used the notation z p = 2NEz p and ( p = 2NE( P . We shall also 
need this kernel for the second argument purely imaginary, 



C(u, i/i) 



uJi(u)I (n) + nJ (u)h(ij) 



v? + jJ? 



(21) 



and in the special case U\ — u 2 — u, 



C{u,u) = - J%(u) + Jf(u) 



f22) 



The major difficulty now is to perform the limit {«/} — > oo in (19). For this 
purpose, it is convenient to divide both numerator and denominator of the 
last term in (19) by A^ f (a) and to perform the limits separately in numerator 
and denominator. Since the derivation is somewhat technical it is presented 
in App. B. The final result is 



where 



and 



det[C(z p ,( q )] P . 



q=l,...,k+N f 



lim 

{a f }^oo det[C(a f , ifig)]f,g=l,...,N f 



det[A 



pq\p,q=l,...,k+Nf 



g=l,...,N f 



(23) 



B = 



-fi>iii(ti>i) 



—fJ>NfIl(fl>N f ) 



— A*l) ' ^-1(^1) •• • (-fJ-N f ) f lNf-l(p>N f 



(24) 



A 



u 



C(ui,ui) 

C(u fc ,ui) 

J (ui) 
UiJi(wi) 

N f- 1 T I \ N f~ 1 T 

1 JNf-l(Ui) ■■■U k J JN f -l{Uk) 



C(ui,u k ) C(ui,i[Xi) ■ ■ -Ciu^ifXNf) 

C(u k ,u k ) C(u k ,ifj,i)---C(u k ,ifji Nf ) 

Jo{u k ) 
u k Ji(u k ) 



B 



(25) 



8 



In compact notation, we have 



A 



pq 



and 



C(u p , u q ) for 1 < p, q < k; 

C(u p , ifi q -k) for 1 < p < k; k + 1 < q < k + Nf] 

uP-^Jp-k-^Uq) for k + 1 < p < k + Nf, 1 < q < k; 

k (-^_ fc )P- fc - 1 /p_ fc _i(^_ fc ) for A; + 1 < p, q < k + N f 

(26) 



Bf 9 = (-ligY^If-iM forl<f,g<N f . (27) 



Rescaling the Rk by (2NE) k , we arrive at the final result for the microscopic 
spectral correlations, 

p k (u^ ...,u k )= ^ R k (x!, ...,x k ) 



k 



I "l r I det[Ap q \p tq= i t ...,k+N f , 

It remains to be shown that S can be identified with the absolute value of the 
chiral condensate (ipip)- According to the Banks-Casher relation [37], we have 
V\ (ipip) | = 7ri?!(0), where i?i(0) is the spectral density at zero. The space-time 
volume V can be identified with 2N. Furthermore, we have 

i?i(0) = 2NE lim pAu) . (29) 

To compute this limit, we use the matrix A in (25) with k = 1 and u± = u. 
We first observe that C(u,u) — > l/(nu) as w — > oo, see (22). In addition, the 
denominator of the entries C(u,ipf), see (21), can be written as a geometric 
series in 1/u 2 so that C(u,ipf) is given as an expansion in 1/u 2 with coeffi- 
cients containing Io(pf), fifli(fif), etc. By subtracting appropriate multiples 
of rows 2 to Nf + 1, the first Nf of these terms can be eliminated, see also 
the discussion of Eq. (B.8) in App. B. The leading large-w behavior of C(u, u) 
is not modified by these subtractions. Higher-order terms in the expansion 
of C(u,ipf) are suppressed by powers of 1/u, the leading term being of or- 
der J(N f +i)mod2(u) /u N f +1 . Even when multiplied by the largest entry in the 
lower-left corner of A, u Nf ~ 1 JN f -i(u), the result is suppressed compared to 
l/{ixu). In the large-w limit, the determinant of (25) with k = 1 thus becomes 
l/(7ru) • det B, and we obtain 

2NE lim pAu) = 2NEu— = -2NE (30) 

and, hence, 2NE = nRi(0) as desired. Therefore, we have shown that the 
functional form of the microscopic spectral correlations in the presence of 
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massive dynamical quarks does not change if a deterministic matrix Y is 
added to the matrix of the Dirac operator, provided that S, i.e., the chiral 
condensate, is nonzero. The only dependence of the final result on the matrix 

Y appears in form of a rescaling of the energy scale, from E at Y = to S at 

Y ^ 0. 

Our final result (28) is given in terms of the determinant of a (k+Nf) x (k+Nf) 
matrix whose entries are simple functions. This structure arises naturally in 
the graded eigenvalue method. Two other forms for the microscopic spec- 
tral correlations have been obtain previously. In the orthogonal-polynomial 
method, the result is given as the determinant of a k x k matrix whose entries 
are (Nf + 2) x (Nf + 2) matrices [17,18]. From the finite volume partition 
function, the result is given as the determinant of a (2k + Nf) x (2k + Nf) 
matrix whose entries are simple functions [20]. At the present time, we do 
not have a closed mathematical proof that these three results are identical. 
However, we have perfomed extensive checks for a large number of different 
values of k and Nf, both numerically and using computer algebra. In all cases, 
the three results agree perfectly so that we do not have any doubt that they 
are identical. One of the virtues of the present method is that it allows us to 
include the deterministic matrix Y as well. Moreover, it appears to lead to the 
most economical representation of the final result. 

In Ref. [18], the distribution of the smallest eigenvalue, P(A m i n ), was also 
computed. Its universality with respect to deformations of Pq(W) was shown in 
Ref. [38]. Since P(A min ) follows directly from the microscopic spectral density 
Pi(u) [18], and since we have shown in this paper that the functional form 
of the latter quantity is not affected by the addition of Y , it follows that the 
functional form of P(A m i n ) also remains unchanged. 



5 Summary 

In this paper, we have extended the graded eigenvalue method to the case 
where massive dynamical quarks enter the distribution of the random matrix. 
The virtue of this approach is twofold. First, we have obtained a novel repre- 
sentation of results computed previously with other methods. Our representa- 
tion appears to be the most economical one. It is also very stable numerically, 
compared to the representations obtained from the orthogonal-polynomial 
method and the finite- volume partition function. Second, our method allows us 
to perform the calculation with a deterministic matrix added to the Dirac ma- 
trix. This is not easily possible in the standard orthogonal-polynomial method 
because this approach rests on the rotation invariance of the matrix ensemble 
which is fully broken due to the presence of the deterministic offset. 
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We point out again that the microscopic correlation functions (28) computed 
in RMT are universal in the sense that they are expected to agree with the 
microscopic spectral correlations of the Dirac operator in full QCD. We hope 
to be able to compare the random matrix results with data from lattice gauge 
simulations with dynamical fermions in the near future. 
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A Derivation of Eq. (7) 



To show that Eq. (7) holds, we write the matrix W in spherical coordinates, 
W = UAV with U G U(iV), V E U(A)/U W (1), and A = diag(Ai, . . . , \ N ), 
where the \ n are real and nonnegative [36]. The integration measure trans- 
forms according to 



d[W] = J(A)d[A]dn(U)dn(V) 



(A.l) 



with 



N 



N 



d[A] = J] d\ n and J(A) = A^(A 2 ) ]J K 



(A.2) 



n=l 



n=l 



The integrations over the \ n extend from to oo, dfj,(U) and dfi{V) are 
the invariant Haar measures, and the Vandermonde determinant is defined 
in Eq. (13). In the integration over d[W] in Eq. (6), we shift W and by 
— Y and obtain 



Z k (J) = 1 f d[A}df,(U)d^V)J(A)e- N ^^ w - Y ^ w '-^ I] det 



irrif —W 
—W^ inif 



x Yl det 

P =i 



x p -W 



x »■ 



Im 



det 



x+ - J p -W 

-w* x;-j p 



-iVE 2 trY 2 



J d[A]F N (A)G N (A)H N (A) 



(A.3) 
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with 



N N f 

F N (A) = A^A 2 ) n Ke- N ^ x ' II ("™/ ~ ^) > ( A -4) 

n=l f=l 

G N (A) = ft ft g " X }_ , (A.5) 

p=l n=l \ X P J P> A n 

H N (A) = J dfi(U)dfi(V)e Nj:2 ^ w+w ^ ¥ . (A.6) 
Analogously, we obtain for Eq. (8) 

-NT, 2 trY 2 . N f N -. 

Z k (J) = 1 -— / d[A]F N (A)G N (A)H N (A) Im II (n+v X 2 ■ ( A - 7 ) 

The normalization factors M and A/" follow by setting (A) to unity in 
Eqs. (A. 3) and (A. 7), respectively. Throughout the remainder of this section, 
we assume that the dummy variables aj are pairwise different. 

We proceed by converting the product over n in (A. 7) to a sum, 
where the subscript (n) means that A n is omitted in A. We now have 



1 _ 71 

(4)^A| = '2a~ f 



Im 7Z+V2 ^ = t 5 ( A ™ ~ °/) + 5 ( A ™ + a /)] 



where in the last step we have used the fact that the integrations over the A n 
extend from to oo and that a/ > (since we consider the limit a/ — > oo). 
Thus, out of the N integrations over the A n in Eq. (A. 7), Nf can be done 
using the (^-functions. Using the symmetry of the integrand with respect to 
the labeling of the A n , we can choose to integrate over the last Nf variables, 
Ajv-iVf+i, • • • , Ajv, by relabeling the A n appropriately and multiplying by a 
combinatorial factor. (The functions -Fat (A) and Gat (A) are obviously sym- 
metric under interchanges A n <-> X m , for Hn(A) this follows from Eq. (A. 14) 
below.) 

We briefly pause to explain the essence of the proof of Eq. (7). After performing 
the Nf integrations in Eq. (A. 7) and taking the limit {a/} — > oo, one obtains 
a result which is essentially equal to the expression for Z k (J) in Eq. (A. 3), the 
only difference being that the integration is over N — Nf variables A n instead of 
over N variables. In the limit N — > oo, this difference can be neglected. (There 
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are some additional prefactors which will be canceled by identical factors in 
the normalization A/".) 

By performing the Nf integrations in Eq. (A. 7) using the (^-functions of Eq. 
(A. 9) the variables Ajv-jv^+i, • • • , Ajv are replaced by a±, . . . , a,N f , respectively. 
We define A' = diag(Ai, . . . , \N-N f ) and A = diag(Ai, . . . , Ajv-jvp oi, • • • , a-N f )- 
The various contributions in the integrand of (A. 7) become 

d[A] — >d[A'}, (A. 10) 

Nf N—Nf N f 

F N (A) — F N - N/ (A')A* N/ (a*) U a f e~ N ^ ft (a} - \ 2 n ) 2 J] (-mj, - a}) 

f=l n=l f'=l 

— > C±(a, m)F N _ Nf (A') as {a/} — > oo , (A. 11) 

fc ^-JV/ T 2 _ \2 A7 2 _ 2 

^(A)^nim n < + p n 7 X 2 n ? + n" 7 2 

p=l n=l - J f) A " /=1 (4 - J p) °/ 

— ► Gj V _jv / (A') as {a/} -> oo . (A. 12) 

Here, Ci(a, to) is a function which no longer depends on the A n . From the 
product over the imaginary parts in Eq. (A. 7) we obtain with Eqs. (A. 8) and 
(A. 9) and after appropriate relabeling of the A n 



f 

Jo 



Nf N 1 



d\ N _ Nf+1 ■ ■ ■ d\ N n im n 



/=i „=i i°/ ) 2 - K 



1 (-7T/2)"/ A! ^ 



— > C 2 (a) as {a/} — > oo , (A. 13) 

where C 2 (a) is a function which depends only on the a/. 

Consider now the angular integrals in Eq. (A. 6). Using Eq. (2.3) of Ref. [36], 
we have 

^ (A) = C MW) ' (A ' 14) 

where we have defined y m = 2NY?y m . The constant c depends neither on A 
nor on Y. After the integration over X N _ Nf+1 , . . . , X N , A is replaced by A, 
i.e., the last Nf entries are ai, . . . , a^ f . We now expand the determinant in 
the numerator of Eq. (A. 14) with respect to the last row. Without loss of 
generality we can assume that y\ < y-i < ... < y^. Using the asymptotic 
behavior of 1$, only the term proportional to Io{a^ f y^) remains in the limit 
a Nf — > oo. The other terms are suppressed by factors of exp(— a^ f (y^ — y n )), 
where n — 1, . . . , N — 1. Using the same argument for the remaining Nf — 1 
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bottom rows, we obtain for {a/} — > oo 

detf/oiAjj^jjn,^!,...^ > d&[Io(^nym)]n,rn=l,...,N-Nf 

x det[/ (a/!/iv-Jv / +a)]/ )9 =i,..,Jv / • (A.15) 
Denoting Y' = diag(yi, . . . , yN-N f ) and Y" = diag(y JV -Jv / +i, . . . , yjv), we have 

A 7 / AT-ATy 

A,v(A 2 ) = A,„ N/ (A' 2 )A N/ (a 2 ) [J II («/ " A «) > (A.16) 

/=1 n=l 

A 7 / AT-ATy 

a^f 2 ) = a^i^a^y" 2 ) n II (y 2 N-N f+f - vl) (A.17) 

f=l n=l 

and thus obtain 

H N (A) — > C 3 (a, Y)H N - Nf (A') as {a f } -> oo . (A.18) 

The function £3(0, Y) no longer depends on the X n . Note that the last Nf 
entries of the diagonal matrix Y have effectively disappeared from the problem. 
However, this effect is negligible in the limit N — > 00. 

Collecting the various terms, we finally obtain for {a/} — > 00 

x I d[A']F JV _ iV/ (A / )G' i v_ J v / (A')//jv_ i v / (A / ) . (A.19) 
The normalization factor A/" is obtained by setting Gjy_]y f (A') to unity so that 



{a/}— >oo 



FN-N f (A')GN-Nf (A')H N - Nf (A') 

lim Z fc (J) = ^ ? (A.20) 

y d[A']F N - Nf (A')H N _ Nf (A') 



which is equal to Z k (J) in Eq. (A. 3) with N replaced by N — Nf. In the limit 
N — > 00, this difference is negligible. This completes the proof. 



B Derivation of Eq. (23) 



We wish to compute the {«/} — > 00 limit in Eq. (19). It is convenient to divide 
both numerator and denominator by the Vandermonde determinant A^ f (a), 
see Eq. (13), and to compute the {«/} — > 00 limit of the quantity 

= det[C(z p , Cg)]p,g=i,...,fc+Ar / ,g ^ 
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with a kernel C given in Eq. (20), z p and ( q defined after Eq. (20), and z p and 
Cj given in Eqs. (15) and (16), respectively. The denominator in Eq. (19) then 
follows immediately by setting k — 0. The a-dependence in the numerator 
determinant of 1Z is found in the last Nf rows, 



1Z = 



A Nf (a) 



C(a 1 ,u 1 ) ■■■ C(a 1 ,u k ) C(a 1 ,ifi 1 ) ■■■ C{ai,iji Nf ) 



C(a Nf ,u 1 ) ■■■ C(a Nf ,u k ) C(a Nf ,i/j 11 ) ■■■ C(a Nf ,ifi Nf ) 



or, indicating rows and their cvdependence schematically by r(«j), 



1Z = 



A Nf (a) 



r(a 1 ) 



r(a Nf 



(B.2) 



Before taking the {«/} — > oo limit, we need to take the limits ctf — > a g for 
all / < g. Let us start with a±. We subtract the row r(a\) from all following 
rows, r(a 2 ) to r(a A r / ), without changing the value of the determinant. The 
Vandermonde determinant supplies factors of the kind 1 /(ccj — a±) for 2 < i < 
Nf. Thus, for the i-th row r(aj), we arrive at 

r(oi) - r(oi) - r( ai ) - ^ ~ ^ ^r(a,) . (B.3) 



Next, the resulting second row <9 Q2 r(a: 2 ) is subtracted from all following rows, 
9a 3 r ( a 3) to d aN ^r(aN f ). Including factors supplied by the Vandermonde de- 
terminant, the j-th row (3 < j < Nf) thus becomes 



d a r(aj) -> d a r(a j ) - d a2 r(a 2 ) 



«2 



(B.4) 



Analogous steps for the subsequent rows lead to higher derivatives, at the 
same time eating up all factors contained in the Vandermonde determinant. 



15 



Finally, we obtain 



lim 1Z = lim 



r(a) 
d a r(a) 



. . . o a r(a) . . . 



(B.5) 



where a is now a simple number. 

We now consider a generic derivative appearing in (B.5). Recall that r(pt) 
stands for C(a,x), where x can be one of the u p or one of the i/if. The 
following manipulations will become more transparent by considering the first 
few special cases. We have 



d° a C(a,x) =C(a,x) 



-> -Ji(a)J (x) 

a 



(B.6) 



The point is that this result is proportional to Jo(x) in the large-a limit. The 
first derivative becomes 



d a C(a,x) 



1 2 

-J (a) ~Ji{ot) 

a a A 

1 



a" 



J Q (x) + —Joi^xJ^x) 



a" 



r J (a)a;Ji(a;) , 



(B.7) 



where the term proportional to Jq(x) has been eliminated by subtracting 
an appropriate multiple of the row (B.6). Thus, the row containing the first 
derivative is proportional to xJ\{x) in the large-a limit. For the second deriva- 
tive, we proceed analogously and obtain, in the large-a limit, terms propor- 
tional to Jo(x), xJi(x), and x 2 J (x). The first two of these terms can be 
eliminated by subtracting appropriate multiples of the rows (B.6) and (B.7), 
respectively. Thus, the row containing the second derivative is proportional to 
x 2 J (x) in the large-a limit. We now proceed to a general m— fold derivative 
for which it is useful to expand the denominator of C in a geometric series, 
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leading to 



d™C(a,x) =£] (™^d l a (-^—^ ST 1 [<*Ji(<*)J (x) - xJoMMx)] 

m fm\ ( 00 r 2n ~ 2 \ 
= E [ l ) d « (E -^rj ST* [«■/!(«) Jb(x) - xJo(a) Mx)] 



m / m \ 00 



Z (2n + Z - l)\x 2n - 2 



S W J S ( ^ (2n-l)! a*"* 



x 



[(aj[ m ' l \a) + (m - /) J^"'" 1 ^)) J (x) - xjt~ l) (a) Ji(x) 

(B.8) 

It will not be necessary to perform the derivatives of Jo (a) and J\{a). The ex- 
pression (B.8) contains terms proportional to Jo(x), xJi(x), ... up to x m J (x) 
(if m is even) or x m J\(x) (if m is odd). All higher order terms in x are 
suppressed by powers of a in the large-a limit, see the sum over n. Sub- 
tracting appropriate multiples of previous rows, only the term proportional to 
x m Jmmod2(x) remains. By using Bessel function recursion relations and adding 
appropriate multiples of previous rows, J mmo d2( x ) can be replaced by J m {x). 
We thus obtain 

d™C(a, x) ^> x m J m {x) x f m (a) (B.9) 
with unspecified functions f m (a). Defining 

N f -1 

•F= Hm J] fm(a) (B.10) 

a °° m=0 

and noting that (ifi) m J m (ifi) = (— fi) m I m (fi), we arrive at 

hm det [^^)U=i,., fc+j v / = ^ detA 
{«/}-«> A Nf (a) 

with A given in Eq. (25). Setting k = 0, we obtain 

hm det[C( ^^ )] (^-^=^-detg (B.12) 

{a/}->oo A A r / (tt) 

with £> given in Eq. (24). Taking the ratio of (B.ll) and (B.12), we finally 
arrive at (23). Note that it is not necessary to evaluate (B.10) since T drops 
out of the final result. 
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